Explaining Biodiversity

Biologists observe wondrous diversity in taxa and traits

Phenotypic Diversity

Morphometricians wish to explain the evolution of phenotypic diversity

Our concern: how do we characterize patterns, and hypothesize processes?

Comparative Biology Tradition

Trait correlations often used to study coevolution and adaptation

Species values commonly utilized

The problem?

Species are not independent of one another

Phylogenetic Comparative Methods (PCMs)

Phylogenetic comparative methods condition the data on the phylogeny during the analysis

Allows one to assess trait covariation while accounting for the non-independence due to shared evolutionary history

PCMs: An Incomplete Historical Walk

The conceptual development of PCMs

Outline

• Phylogenetic regression/anova/association models

  • Phylogenetic Generalized Least Squares methods

  • Phylogenetic PLS

• Phylogenetic ordination

  • Phylomorphospace

• Exploring and modelling evolutionary processes

  • Phylogenetic signal

  • Evolutionary rates

  • Evolutionary models

PCMs: General Statistical Concepts

• PCMs condition the data on the phylogeny under an evolutionary model

• Data conditioning: Account for phylogenetic non-independence

  • (i.e., examine pattern ‘given’ the phylogeny)

• Evolutionary model: How trait variation is expected to accumulate

  • Brownian motion (BM): A common null model of trait change
  • No change in \(\small\mu\) but \(\small\sigma^2\uparrow\propto{time}\)

PCM: General Model

Most PCMs use GLS (generalized least squares) as a model:

\[\small\mathbf{Y}=\mathbf{X{\hat{\beta}}+\epsilon}\] Here, \(\small\epsilon\) is not iid but is \(\small\sim\mathcal{N}(0,\textbf{V})\): containing expected covariation between taxa as described by the phylogeny (in matrix \(\small\mathbf{V}\))

\(\small\mathbf{V}\) is the phylogenetic covariance matrix

Describes the amount of evolutionary time species share via common ancestry (and thus how similar they are expected to be)

Sometimes V is called C (particularly when derived under Brownian motion)

PCM: General Model

\[\small\mathbf{Y}=\mathbf{X{\hat{\beta}}+\epsilon}\]

Model design (\(\small\mathbf{X}\)) describes the type of analysis

Parameters of model (and model significance) obtained in various ways

see Adams and Collyer. Syst. Biol. (2018); Adams and Collyer. Ann. Rev. Ecol. Evol. Syst. (2019)

1: Phylogenetic Regression

Evaluate \(\small\mathbf{Y}=\mathbf{X}\mathbf{\beta } +\mathbf{E}\) in a phylogenetic context

The workhorse of PCMs

1: Phylogenetic Regression

ANOVA and regression models that account for phylogeny

Requires a model of evolutionary change: typically Brownian motion (BM)

1: Phylogenetic Regression

ANOVA and regression models that account for phylogeny

Requires a model of evolutionary change: typically Brownian motion (BM)

\(\Delta\mu = 0\)

\(\sigma^2\) (variance among taxa) \(\uparrow\) \(\propto\) time

1: Phylogenetic Regression

ANOVA and regression models that account for phylogeny

Requires a model of evolutionary change: typically Brownian motion (BM)

\(\Delta\mu = 0\)

\(\sigma^2\) (variance among taxa) \(\uparrow\) \(\propto\) time

Statistical model: \(\small\mathbf{Y}=\mathbf{X}\mathbf{\beta } +\mathbf{E}\) where \(\small\mathbf{E} \sim \mathcal{N}(0,\textbf{V})\)

Several implementations possible (yield identical results if implemented properly)

Phylogenetically Independent Contrasts

Estimate contrast scores between pairs of taxa (tips or nodes)

Use contrasts for analyses (OLS solution)

see Felsenstein. Am. Nat. (1985)

Phylogenetically Independent Contrasts

Coefficients found as: \(\small\hat{\mathbf{\beta }}=\left ( \mathbf{X}^{T}_{pic} \mathbf{X}_{pic}\right )^{-1}\left ( \mathbf{X}^{T}_{pic} \mathbf{Y}_{pic}\right )\)

Phylogenetically Independent Contrasts: Example

What is the association between Y and X?

Phylogenetically Independent Contrasts: Example

## Registered S3 method overwritten by 'geiger':
##   method            from
##   unique.multiPhylo ape

What is the association between Y and X?

## Analysis of Variance Table
## 
## Response: pic.y
##           Df  Sum Sq Mean Sq F value  Pr(>F)   
## pic.x      1 14.3519 14.3519  19.285 0.00461 **
## Residuals  6  4.4651  0.7442                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

##     pic.x 
## 0.8846971

Phylogenetic Generalized Least Squares (PGLS)

Use GLS model with non-independent error structure

Statistical model: \(\small\mathbf{Y}=\mathbf{X}\mathbf{\beta } +\mathbf{E}\) where \(\small\mathbf{E} \sim \mathcal{N}(0,\textbf{V})\)

Phylogenetic Generalized Least Squares (PGLS)

Use GLS model with non-independent error structure

Statistical model: \(\small\mathbf{Y}=\mathbf{X}\mathbf{\beta } +\mathbf{E}\) where \(\small\mathbf{E} \sim \mathcal{N}(0,\textbf{V})\)

\(\textbf{V}\) is the phylogenetic covariance matrix

Describes expected covariance among taxa due to shared evolutionary history (typically under BM)

Coefficients found as: \(\small\hat{\mathbf{\beta }}=\left ( \mathbf{X}^{T} \mathbf{V}^{-1} \mathbf{X}\right )^{-1}\left ( \mathbf{X}^{T} \mathbf{V}^{-1}\mathbf{Y}\right )\)

see Grafen. Phil. Trans. Roy. Soc. (1989)

PGLS: Example

## Denom. DF: 6 
##             numDF  F-value p-value
## (Intercept)     1 12.87792  0.0115
## X               1 19.28544  0.0046

## (Intercept)           X 
##  -2.3296557   0.8846971

Identical results to PICs!

Statistical Digression: OLS vs GLS

To understand what PGLS is doing, consider the following methods for obtaining model parameters

OLS model: \(\small\hat{\mathbf{\beta }}=\left ( \mathbf{X}^{T}\mathbf{X}\right )^{-1}\left ( \mathbf{X}^{T} \mathbf{Y}\right )\)

Statistical Digression: OLS vs GLS

To understand what PGLS is doing, consider the following methods for obtaining model parameters

OLS model: \(\small\hat{\mathbf{\beta }}=\left ( \mathbf{X}^{T}\mathbf{X}\right )^{-1}\left ( \mathbf{X}^{T} \mathbf{Y}\right )\)

Unweighted GLS model: \(\small\hat{\mathbf{\beta }}=\left ( \mathbf{X}^{T} \mathbf{V}_{id}^{-1} \mathbf{X}\right )^{-1}\left ( \mathbf{X}^{T} \mathbf{V}_{id}^{-1}\mathbf{Y}\right )\) where

\[\tiny\mathbf{V}_{id}= \left( \begin{array}{ccc} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{array} \right) \]

Statistical Digression: OLS vs GLS

To understand what PGLS is doing, consider the following methods for obtaining model parameters

OLS model: \(\small\hat{\mathbf{\beta }}=\left ( \mathbf{X}^{T}\mathbf{X}\right )^{-1}\left ( \mathbf{X}^{T} \mathbf{Y}\right )\)

Unweighted GLS model: \(\small\hat{\mathbf{\beta }}=\left ( \mathbf{X}^{T} \mathbf{V}_{id}^{-1} \mathbf{X}\right )^{-1}\left ( \mathbf{X}^{T} \mathbf{V}_{id}^{-1}\mathbf{Y}\right )\) where

\[\tiny\mathbf{V}_{id}= \left( \begin{array}{ccc} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{array} \right) \]

Weighted GLS model: \(\small\hat{\mathbf{\beta }}=\left ( \mathbf{X}^{T} \mathbf{V}_{phy}^{-1} \mathbf{X}\right )^{-1}\left ( \mathbf{X}^{T} \mathbf{V}_{phy}^{-1}\mathbf{Y}\right )\) where

\[\tiny\mathbf{V}_{id}= \left( \begin{array}{ccc} (v_1+v_2) & v_{12} & 0 \\ v_{12} & (v_1+v_2) & 0 \\ 0 & 0 & 1 \end{array} \right) \]

In PGLS, the weights are the phylogenetic distances, which describe the phylogenetic non-independence

Statistical Digression: OLS vs GLS

To understand what PGLS is doing, consider the following methods for obtaining model parameters

OLS model: \(\small\hat{\mathbf{\beta }}=\left ( \mathbf{X}^{T}\mathbf{X}\right )^{-1}\left ( \mathbf{X}^{T} \mathbf{Y}\right )\)

Unweighted GLS model: \(\small\hat{\mathbf{\beta }}=\left ( \mathbf{X}^{T} \mathbf{V}_{id}^{-1} \mathbf{X}\right )^{-1}\left ( \mathbf{X}^{T} \mathbf{V}_{id}^{-1}\mathbf{Y}\right )\) where

\[\tiny\mathbf{V}_{id}= \left( \begin{array}{ccc} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{array} \right) \]

Weighted GLS model: \(\small\hat{\mathbf{\beta }}=\left ( \mathbf{X}^{T} \mathbf{V}_{phy}^{-1} \mathbf{X}\right )^{-1}\left ( \mathbf{X}^{T} \mathbf{V}_{phy}^{-1}\mathbf{Y}\right )\) where

\[\tiny\mathbf{V}_{id}= \left( \begin{array}{ccc} (v_1+v_2) & v_{12} & 0 \\ v_{12} & (v_1+v_2) & 0 \\ 0 & 0 & 1 \end{array} \right) \]

In PGLS, the weights are the phylogenetic distances, which describe the phylogenetic non-independence

Phylogenetic Transformation

Utilize OLS transformation of GLS models

Phylogenetic transformation matrix: \(\small\mathbf{T}=\left ( \mathbf{U}\mathbf{W}^{-1/2} \mathbf{U}^{T}\right )^{-1}\)

U and W are eigenvectors and eigenvalues of V

Transformed data: \(\small\tilde{\mathbf{Y}}=\mathbf{TY}\)

Transformed design matrix: \(\small\tilde{\mathbf{X}}=\mathbf{TX}\)

Coefficients solved as: \(\small\hat{\mathbf{\beta }}=\left ( \mathbf{\tilde{X}}^{T} \mathbf{\tilde{X}}\right )^{-1}\left ( \mathbf{\tilde{X}}^{T} \mathbf{\tilde{Y}}\right )\)

see Garland and Ives. Am. Nat. (2000); Adams. Evol. (2014); Adams and Collyer. Evol. (2018)

Phylogenetic Transformation: Example

##           Df      SS      MS     Rsq      F      Z Pr(>F)   
## X          1 14.3519 14.3519 0.76271 19.285 1.8228  0.002 **
## Residuals  6  4.4651  0.7442 0.23729                        
## Total      7 18.8170                                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## (Intercept)           X 
##  -2.3296557   0.8846971

Identical results to PICs & PGLS!

Assessing Significance

PIC, PGLS, and Phylogenetic transform yield identical parameter estimates

\[\tiny \hat{\mathbf{\beta }}=\left ( \mathbf{X}^{T} \mathbf{V}^{-1} \mathbf{X}\right )^{-1}\left ( \mathbf{X}^{T} \mathbf{V}^{-1}\mathbf{Y}\right )\]

\[\tiny =\left ( \mathbf{X}^{T}_{pic} \mathbf{X}_{pic}\right )^{-1}\left ( \mathbf{X}^{T}_{pic} \mathbf{Y}_{pic}\right )\]

\[\tiny =\left ( \mathbf{\tilde{X}}^{T} \mathbf{\tilde{X}}\right )^{-1}\left ( \mathbf{\tilde{X}}^{T} \mathbf{\tilde{Y}}\right )\]

• Model significance typically accomplished using parametric procedures

  • F-ratios compared to parametric F-distribution

  • Optimize \(\log{L}\) for model

\[\tiny \log{L}=\log{ \left(\frac{exp{\left({-\frac{1}{2}{\left({\mathbf{Y}-E(\mathbf{Y})}\right)^{T}}} \mathbf{V}^{-1}\left({\mathbf{Y}-E(\mathbf{Y})}\right)\right)}} {\sqrt{\left({2\pi}\right)^{Np}\times\begin{vmatrix} \mathbf{V} \end{vmatrix}}}\right)}\]

Assessing Significance

PROBLEM: Parametric significance testing suffers from Rao’s paradox

Power reduces as p-dimensions increase

Another solution required

Adams. Evol. (2014); Adams and Collyer. Syst. Biol. (2018)

Assessing Significance: RRPP

Transform data: \(\small\tilde{\mathbf{Y}}=\mathbf{TY}\), \(\small\tilde{\mathbf{X}}_{F}=\mathbf{TX}_{F}\), \(\small\tilde{\mathbf{X}}_{R}=\mathbf{TX}_{R}\)

Obtain Parameter estimates: \(\small\hat{\mathbf{\beta }}=\left ( \mathbf{\tilde{X_{F}}}^{T} \mathbf{\tilde{X_{F}}}\right )^{-1}\left ( \mathbf{\tilde{X_{F}}}^{T} \mathbf{\tilde{Y}}\right )\), and permute residuals (as described previously)

Assessing Significance: RRPP

Transform data: \(\small\tilde{\mathbf{Y}}=\mathbf{TY}\), \(\small\tilde{\mathbf{X}}_{F}=\mathbf{TX}_{F}\), \(\small\tilde{\mathbf{X}}_{R}=\mathbf{TX}_{R}\)

Obtain Parameter estimates: \(\small\hat{\mathbf{\beta }}=\left ( \mathbf{\tilde{X_{F}}}^{T} \mathbf{\tilde{X_{F}}}\right )^{-1}\left ( \mathbf{\tilde{X_{F}}}^{T} \mathbf{\tilde{Y}}\right )\), and permute residuals (as described previously)

Appropriate for ANOVA, regression and more complex PGLS models with high-dimensional data (Adams & Collyer. Evol. 2018)

RULE: Phylo-transform first, permute second!

Assessing Significance: RRPP

Transform data: \(\small\tilde{\mathbf{Y}}=\mathbf{TY}\), \(\small\tilde{\mathbf{X}}_{F}=\mathbf{TX}_{F}\), \(\small\tilde{\mathbf{X}}_{R}=\mathbf{TX}_{R}\)

Obtain Parameter estimates: \(\small\hat{\mathbf{\beta }}=\left ( \mathbf{\tilde{X_{F}}}^{T} \mathbf{\tilde{X_{F}}}\right )^{-1}\left ( \mathbf{\tilde{X_{F}}}^{T} \mathbf{\tilde{Y}}\right )\), and permute residuals (as described previously)

Appropriate for ANOVA, regression and more complex PGLS models with high-dimensional data (Adams & Collyer. Evol. 2018)

RULE: Phylo-transform first, permute second!

Note: alternatives permuting PICs (Klingenberg and and Marugan-Loban 2013) are not correct (see Adams and Collyer 2015. Evolution).

WHAT YOU SHUFFLE MATTERS! (must shuffle correct exchangeable units: see Adams and Collyer. Evol. 2015; Adams and Collyer. Evol. 2018)

Phylogenetic Regression: Example

Does head shape covary with body size among Plethodon salamander species?

Phylogenetic Regression: Example

Are body dimensions associated with overall size across Plethodon salamander species?

##           Df        SS         MS     Rsq      F      Z Pr(>Cohen's f-squared)
## svl        1 0.0006586 0.00065861 0.07039 3.0289 2.0186                  0.012
## Residuals 40 0.0086977 0.00021744 0.92961                                     
## Total     41 0.0093563                                                        
##            
## svl       *
## Residuals  
## Total      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

YES! There is a significant association between head shape and body size in salamanders

Phylogenetic ANOVA: Example

Does head shape differ between high-elevation and low-elevation salamander species?

Phylogenetic ANOVA: Example

Does head shape differ between high-elevation and low-elevation salamander species?

##           Df        SS         MS     Rsq      F      Z Pr(>Cohen's f-squared)
## elev       1 0.0006142 0.00061420 0.06565 2.8103 1.8621                  0.033
## Residuals 40 0.0087421 0.00021855 0.93435                                     
## Total     41 0.0093563                                                        
##            
## elev      *
## Residuals  
## Total      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

YES! Head shape differs in high- and low-elevation species

Phylogenetic ANOVA: Visualizations

One can visualize group dispersion in morphospace…

BROKEN: don’t know why

#{r echo = FALSE, eval=TRUE,out.width="80%"} #plot.res <- gm.prcomp(shape,phy=plethtree, data=gdf) #ord.plot <- plot(plot.res,phylo = FALSE, pch=21, bg=gdf$elev, cex=1.5) #shapeHulls(ord.plot, groups = gdf$elev, # group.cols = c("red", "black"), # group.lwd = rep(1, 2), group.lty = c(2, 1)) #legend("topright", levels(gdf$elev), # col = c("black", "red"), # lwd = rep(1,2), lty = c(2, 1)) #

Phylogenetic ANOVA: Visualizations

… and via TPS deformation grids

Phylogenetic ANOVA: Group Aggregation

How groups are distributed on the phylogeny can make a difference: Adams and Collyer. Evol. (2018)

Too few group ‘state’ changes across phylogeny results in lower power